module FrictionVelocityMod

!------------------------------------------------------------------------------
!BOP
!
! !MODULE: FrictionVelocityMod
!
! !DESCRIPTION:
! Calculation of the friction velocity, relation for potential
! temperature and humidity profiles of surface boundary layer.
!
! !USES:
  use shr_kind_mod, only: r8 => shr_kind_r8
!
! !PUBLIC TYPES:
  implicit none
  save
!
! !PUBLIC MEMBER FUNCTIONS:
  public :: FrictionVelocity       ! Calculate friction velocity
  public :: MoninObukIni           ! Initialization of the Monin-Obukhov length
!
! !PRIVATE MEMBER FUNCTIONS:
  private :: StabilityFunc1        ! Stability function for rib < 0.
  private :: StabilityFunc2        ! Stability function for rib < 0.
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!------------------------------------------------------------------------------

contains

!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: FrictionVelocity
!
! !INTERFACE:
  subroutine FrictionVelocity(lbn, ubn, fn, filtern, &
                              displa, z0m, z0h, z0q, &
                              obu, iter, ur, um, ustar, &
                              temp1, temp2, temp12m, temp22m, fm, landunit_index)
!
! !DESCRIPTION:
! Calculation of the friction velocity, relation for potential
! temperature and humidity profiles of surface boundary layer.
! The scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
! Vol. 11, 2628-2644.
!
! !USES:
   use clmtype
   use clm_atmlnd, only : clm_a2l
   use clm_varcon, only : vkc
   use clm_varctl, only : iulog
!
! !ARGUMENTS:
   implicit none
   integer , intent(in)  :: lbn, ubn         ! pft/landunit array bounds
   integer , intent(in)  :: fn               ! number of filtered pft/landunit elements
   integer , intent(in)  :: filtern(fn)      ! pft/landunit filter
   real(r8), intent(in)  :: displa(lbn:ubn)  ! displacement height (m)
   real(r8), intent(in)  :: z0m(lbn:ubn)     ! roughness length over vegetation, momentum [m]
   real(r8), intent(in)  :: z0h(lbn:ubn)     ! roughness length over vegetation, sensible heat [m]
   real(r8), intent(in)  :: z0q(lbn:ubn)     ! roughness length over vegetation, latent heat [m]
   real(r8), intent(in)  :: obu(lbn:ubn)     ! monin-obukhov length (m)
   integer,  intent(in)  :: iter             ! iteration number
   real(r8), intent(in)  :: ur(lbn:ubn)      ! wind speed at reference height [m/s]
   real(r8), intent(in)  :: um(lbn:ubn)      ! wind speed including the stablity effect [m/s]
   logical,  optional, intent(in)  :: landunit_index  ! optional argument that defines landunit or pft level
   real(r8), intent(out) :: ustar(lbn:ubn)   ! friction velocity [m/s]
   real(r8), intent(out) :: temp1(lbn:ubn)   ! relation for potential temperature profile
   real(r8), intent(out) :: temp12m(lbn:ubn) ! relation for potential temperature profile applied at 2-m
   real(r8), intent(out) :: temp2(lbn:ubn)   ! relation for specific humidity profile
   real(r8), intent(out) :: temp22m(lbn:ubn) ! relation for specific humidity profile applied at 2-m
   real(r8), intent(inout) :: fm(lbn:ubn)    ! diagnose 10m wind (DUST only)
!
! !CALLED FROM:
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
! 12/19/01, Peter Thornton
! Added arguments to eliminate passing clm derived type into this function.
! Created by Mariana Vertenstein
!
! !LOCAL VARIABLES:
!
! local pointers to implicit in arguments
!
   integer , pointer :: ngridcell(:)      !pft/landunit gridcell index
   real(r8), pointer :: forc_hgt_u_pft(:) !observational height of wind at pft level [m]
   real(r8), pointer :: forc_hgt_t_pft(:) !observational height of temperature at pft level [m]
   real(r8), pointer :: forc_hgt_q_pft(:) !observational height of specific humidity at pft level [m]
   integer , pointer :: pfti(:)           !beginning pfti index for landunit
   integer , pointer :: pftf(:)           !final pft index for landunit
!
! local pointers to implicit out arguments
!
   real(r8), pointer :: u10(:)         ! 10-m wind (m/s) (for dust model)
   real(r8), pointer :: fv(:)          ! friction velocity (m/s) (for dust model)
   real(r8), pointer :: vds(:)         ! dry deposition velocity term (m/s) (for SO4 NH4NO3)
   real(r8), pointer :: u10_clm(:)     ! 10-m wind (m/s)
   real(r8), pointer :: va(:)          ! atmospheric wind speed plus convective velocity (m/s)
!
!
! !OTHER LOCAL VARIABLES:
!EOP
!
   real(r8), parameter :: zetam = 1.574_r8 ! transition point of flux-gradient relation (wind profile)
   real(r8), parameter :: zetat = 0.465_r8 ! transition point of flux-gradient relation (temp. profile)
   integer :: f                            ! pft/landunit filter index
   integer :: n                            ! pft/landunit index
   integer :: g                            ! gridcell index
   integer :: pp                           ! pfti,pftf index
   real(r8):: zldis(lbn:ubn)               ! reference height "minus" zero displacement heght [m]
   real(r8):: zeta(lbn:ubn)                ! dimensionless height used in Monin-Obukhov theory
   real(r8) :: tmp1,tmp2,tmp3,tmp4         ! Used to diagnose the 10 meter wind
   real(r8) :: fmnew                       ! Used to diagnose the 10 meter wind
   real(r8) :: fm10                        ! Used to diagnose the 10 meter wind
   real(r8) :: zeta10                      ! Used to diagnose the 10 meter wind
   real(r8) :: vds_tmp                     ! Temporary for dry deposition velocity
!------------------------------------------------------------------------------

   ! Assign local pointers to derived type members (gridcell-level)

   if (present(landunit_index)) then
     ngridcell  => clm3%g%l%gridcell
   else
     ngridcell  => clm3%g%l%c%p%gridcell
   end if

   vds        => clm3%g%l%c%p%pps%vds
   u10        => clm3%g%l%c%p%pps%u10
   u10_clm    => clm3%g%l%c%p%pps%u10_clm
   va         => clm3%g%l%c%p%pps%va
   fv         => clm3%g%l%c%p%pps%fv

   ! Assign local pointers to derived type members (pft or landunit-level)

   pfti             => clm3%g%l%pfti
   pftf             => clm3%g%l%pftf

   ! Assign local pointers to derived type members (pft-level)

   forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft
   forc_hgt_t_pft => clm3%g%l%c%p%pps%forc_hgt_t_pft
   forc_hgt_q_pft => clm3%g%l%c%p%pps%forc_hgt_q_pft

   ! Adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions.

#if (!defined PERGRO)

   do f = 1, fn
      n = filtern(f)
      g = ngridcell(n)

      ! Wind profile

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_u_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_u_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetam) then
         ustar(n) = vkc*um(n)/(log(-zetam*obu(n)/z0m(n))&
              - StabilityFunc1(-zetam) &
              + StabilityFunc1(z0m(n)/obu(n)) &
              + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8))
      else if (zeta(n) < 0._r8) then
         ustar(n) = vkc*um(n)/(log(zldis(n)/z0m(n))&
              - StabilityFunc1(zeta(n))&
              + StabilityFunc1(z0m(n)/obu(n)))
      else if (zeta(n) <=  1._r8) then
         ustar(n) = vkc*um(n)/(log(zldis(n)/z0m(n)) + 5._r8*zeta(n) -5._r8*z0m(n)/obu(n))
      else
         ustar(n) = vkc*um(n)/(log(obu(n)/z0m(n))+5._r8-5._r8*z0m(n)/obu(n) &
              +(5._r8*log(zeta(n))+zeta(n)-1._r8))
      end if
      
      if (zeta(n) < 0._r8) then
         vds_tmp = 2.e-3_r8*ustar(n) * ( 1._r8 + (300._r8/(-obu(n)))**0.666_r8)
      else
         vds_tmp = 2.e-3_r8*ustar(n)
      endif

      if (present(landunit_index)) then
         do pp = pfti(n),pftf(n)
            vds(pp) = vds_tmp
         end do
      else
         vds(n) = vds_tmp
      end if

! Calculate a 10-m wind (10m + z0m + d)
! For now, this will not be the same as the 10-m wind calculated for the dust 
! model because the CLM stability functions are used here, not the LSM stability
! functions used in the dust model. We will eventually change the dust model to be 
! consistent with the following formulation.
! Note that the 10-m wind calculated this way could actually be larger than the
! atmospheric forcing wind because 1) this includes the convective velocity, 2)
! this includes the 1 m/s minimum wind threshold

! If forcing height is less than or equal to 10m, then set 10-m wind to um
      if (present(landunit_index)) then
        do pp = pfti(n),pftf(n)
          if (zldis(n)-z0m(n) .le. 10._r8) then
            u10_clm(pp) = um(n)
          else
            if (zeta(n) < -zetam) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(-zetam*obu(n)/(10._r8+z0m(n)))      &
                                      - StabilityFunc1(-zetam)                              &
                                      + StabilityFunc1((10._r8+z0m(n))/obu(n))              &
                                      + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8)) )
            else if (zeta(n) < 0._r8) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n)))           &
                                      - StabilityFunc1(zeta(n))                             &
                                      + StabilityFunc1((10._r8+z0m(n))/obu(n))) )
            else if (zeta(n) <=  1._r8) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n)))           &
                                      + 5._r8*zeta(n) - 5._r8*(10._r8+z0m(n))/obu(n)) )
            else
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(obu(n)/(10._r8+z0m(n)))             &
                                      + 5._r8 - 5._r8*(10._r8+z0m(n))/obu(n)                &
                                      + (5._r8*log(zeta(n))+zeta(n)-1._r8)) )

            end if
          end if
          va(pp) = um(n)
        end do
      else
        if (zldis(n)-z0m(n) .le. 10._r8) then
          u10_clm(n) = um(n)
        else
          if (zeta(n) < -zetam) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(-zetam*obu(n)/(10._r8+z0m(n)))         &
                                   - StabilityFunc1(-zetam)                                 &
                                   + StabilityFunc1((10._r8+z0m(n))/obu(n))                 &
                                   + 1.14_r8*((-zeta(n))**0.333_r8-(zetam)**0.333_r8)) )
          else if (zeta(n) < 0._r8) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n)))              &
                                   - StabilityFunc1(zeta(n))                                &
                                   + StabilityFunc1((10._r8+z0m(n))/obu(n))) )
          else if (zeta(n) <=  1._r8) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n)))              &
                                   + 5._r8*zeta(n) - 5._r8*(10._r8+z0m(n))/obu(n)) )
          else
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(obu(n)/(10._r8+z0m(n)))                &
                                   + 5._r8 - 5._r8*(10._r8+z0m(n))/obu(n)                   &
                                   + (5._r8*log(zeta(n))+zeta(n)-1._r8)) )
          end if
        end if
        va(n) = um(n)
      end if

      ! Temperature profile

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_t_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_t_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp1(n) = vkc/(log(-zetat*obu(n)/z0h(n))&
              - StabilityFunc2(-zetat) &
              + StabilityFunc2(z0h(n)/obu(n)) &
              + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
      else if (zeta(n) < 0._r8) then
         temp1(n) = vkc/(log(zldis(n)/z0h(n)) &
              - StabilityFunc2(zeta(n)) &
              + StabilityFunc2(z0h(n)/obu(n)))
      else if (zeta(n) <=  1._r8) then
         temp1(n) = vkc/(log(zldis(n)/z0h(n)) + 5._r8*zeta(n) - 5._r8*z0h(n)/obu(n))
      else
         temp1(n) = vkc/(log(obu(n)/z0h(n)) + 5._r8 - 5._r8*z0h(n)/obu(n) &
              + (5._r8*log(zeta(n))+zeta(n)-1._r8))
      end if

      ! Humidity profile

      if (present(landunit_index)) then
       if (forc_hgt_q_pft(pfti(n)) == forc_hgt_t_pft(pfti(n)) .and. z0q(n) == z0h(n)) then
         temp2(n) = temp1(n)
       else
         zldis(n) = forc_hgt_q_pft(pfti(n))-displa(n)
         zeta(n) = zldis(n)/obu(n)
         if (zeta(n) < -zetat) then
            temp2(n) = vkc/(log(-zetat*obu(n)/z0q(n)) &
                 - StabilityFunc2(-zetat) &
                 + StabilityFunc2(z0q(n)/obu(n)) &
                 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
         else if (zeta(n) < 0._r8) then
            temp2(n) = vkc/(log(zldis(n)/z0q(n)) &
                 - StabilityFunc2(zeta(n)) &
                 + StabilityFunc2(z0q(n)/obu(n)))
         else if (zeta(n) <=  1._r8) then
            temp2(n) = vkc/(log(zldis(n)/z0q(n)) + 5._r8*zeta(n)-5._r8*z0q(n)/obu(n))
         else
            temp2(n) = vkc/(log(obu(n)/z0q(n)) + 5._r8 - 5._r8*z0q(n)/obu(n) &
                 + (5._r8*log(zeta(n))+zeta(n)-1._r8))
         end if
       end if
      else
       if (forc_hgt_q_pft(n) == forc_hgt_t_pft(n) .and. z0q(n) == z0h(n)) then
         temp2(n) = temp1(n)
       else
         zldis(n) = forc_hgt_q_pft(n)-displa(n)
         zeta(n) = zldis(n)/obu(n)
         if (zeta(n) < -zetat) then
            temp2(n) = vkc/(log(-zetat*obu(n)/z0q(n)) &
                 - StabilityFunc2(-zetat) &
                 + StabilityFunc2(z0q(n)/obu(n)) &
                 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
         else if (zeta(n) < 0._r8) then
            temp2(n) = vkc/(log(zldis(n)/z0q(n)) &
                 - StabilityFunc2(zeta(n)) &
                 + StabilityFunc2(z0q(n)/obu(n)))
         else if (zeta(n) <=  1._r8) then
            temp2(n) = vkc/(log(zldis(n)/z0q(n)) + 5._r8*zeta(n)-5._r8*z0q(n)/obu(n))
         else
            temp2(n) = vkc/(log(obu(n)/z0q(n)) + 5._r8 - 5._r8*z0q(n)/obu(n) &
                 + (5._r8*log(zeta(n))+zeta(n)-1._r8))
         end if
       endif
      endif

      ! Temperature profile applied at 2-m

      zldis(n) = 2.0_r8 + z0h(n)
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp12m(n) = vkc/(log(-zetat*obu(n)/z0h(n))&
              - StabilityFunc2(-zetat) &
              + StabilityFunc2(z0h(n)/obu(n)) &
              + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
      else if (zeta(n) < 0._r8) then
         temp12m(n) = vkc/(log(zldis(n)/z0h(n)) &
              - StabilityFunc2(zeta(n))  &
              + StabilityFunc2(z0h(n)/obu(n)))
      else if (zeta(n) <=  1._r8) then
         temp12m(n) = vkc/(log(zldis(n)/z0h(n)) + 5._r8*zeta(n) - 5._r8*z0h(n)/obu(n))
      else
         temp12m(n) = vkc/(log(obu(n)/z0h(n)) + 5._r8 - 5._r8*z0h(n)/obu(n) &
              + (5._r8*log(zeta(n))+zeta(n)-1._r8))
      end if

      ! Humidity profile applied at 2-m

      if (z0q(n) == z0h(n)) then
         temp22m(n) = temp12m(n)
      else
         zldis(n) = 2.0_r8 + z0q(n)
         zeta(n) = zldis(n)/obu(n)
         if (zeta(n) < -zetat) then
            temp22m(n) = vkc/(log(-zetat*obu(n)/z0q(n)) - &
                 StabilityFunc2(-zetat) + StabilityFunc2(z0q(n)/obu(n)) &
                 + 0.8_r8*((zetat)**(-0.333_r8)-(-zeta(n))**(-0.333_r8)))
         else if (zeta(n) < 0._r8) then
            temp22m(n) = vkc/(log(zldis(n)/z0q(n)) - &
                 StabilityFunc2(zeta(n))+StabilityFunc2(z0q(n)/obu(n)))
         else if (zeta(n) <=  1._r8) then
            temp22m(n) = vkc/(log(zldis(n)/z0q(n)) + 5._r8*zeta(n)-5._r8*z0q(n)/obu(n))
         else
            temp22m(n) = vkc/(log(obu(n)/z0q(n)) + 5._r8 - 5._r8*z0q(n)/obu(n) &
                 + (5._r8*log(zeta(n))+zeta(n)-1._r8))
         end if
      end if

      ! diagnose 10-m wind for dust model (dstmbl.F)
      ! Notes from C. Zender's dst.F:
      ! According to Bon96 p. 62, the displacement height d (here displa) is
      ! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
      ! Therefore d <= 0.034*z1 and may safely be neglected.
      ! Code from LSM routine SurfaceTemperature was used to obtain u10

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_u_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_u_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (min(zeta(n), 1._r8) < 0._r8) then
         tmp1 = (1._r8 - 16._r8*min(zeta(n),1._r8))**0.25_r8
         tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
         tmp3 = log((1._r8+tmp1)/2._r8)
         fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
      else
         fmnew = -5._r8*min(zeta(n),1._r8)
      endif
      if (iter == 1) then
         fm(n) = fmnew
      else
         fm(n) = 0.5_r8 * (fm(n)+fmnew)
      end if
      zeta10 = min(10._r8/obu(n), 1._r8)
      if (zeta(n) == 0._r8) zeta10 = 0._r8
      if (zeta10 < 0._r8) then
         tmp1 = (1.0_r8 - 16.0_r8 * zeta10)**0.25_r8
         tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
         tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
         fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
      else                ! not stable
         fm10 = -5.0_r8 * zeta10
      end if
      if (present(landunit_index)) then
        tmp4 = log( max( 1.0_8, forc_hgt_u_pft(pfti(n)) / 10._r8) )
      else 
        tmp4 = log( max( 1.0_8, forc_hgt_u_pft(n) / 10._r8) )
      end if
      if (present(landunit_index)) then
        do pp = pfti(n),pftf(n)
          u10(pp) = ur(n) - ustar(n)/vkc * (tmp4 - fm(n) + fm10)
          fv(pp)  = ustar(n)
        end do 
      else
        u10(n) = ur(n) - ustar(n)/vkc * (tmp4 - fm(n) + fm10)
        fv(n)  = ustar(n)
      end if

   end do
#endif


#if (defined PERGRO)

   !===============================================================================
   ! The following only applies when PERGRO is defined
   !===============================================================================

   do f = 1, fn
      n = filtern(f)
      g = ngridcell(n)

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_u_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_u_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetam) then           ! zeta < -1
         ustar(n) = vkc * um(n) / log(-zetam*obu(n)/z0m(n))
      else if (zeta(n) < 0._r8) then         ! -1 <= zeta < 0
         ustar(n) = vkc * um(n) / log(zldis(n)/z0m(n))
      else if (zeta(n) <= 1._r8) then        !  0 <= ztea <= 1
         ustar(n)=vkc * um(n)/log(zldis(n)/z0m(n))
      else                             !  1 < zeta, phi=5+zeta
         ustar(n)=vkc * um(n)/log(obu(n)/z0m(n))
      endif

! Calculate a 10-m wind (10m + z0m + d)
! For now, this will not be the same as the 10-m wind calculated for the dust 
! model because the CLM stability functions are used here, not the LSM stability
! functions used in the dust model. We will eventually change the dust model to be 
! consistent with the following formulation.
! Note that the 10-m wind calculated this way could actually be larger than the
! atmospheric forcing wind because 1) this includes the convective velocity, 2)
! this includes the 1 m/s minimum wind threshold

! If forcing height is less than or equal to 10m, then set 10-m wind to um
      if (present(landunit_index)) then
        do pp = pfti(n),pftf(n)
          if (zldis(n)-z0m(n) .le. 10._r8) then
            u10_clm(pp) = um(n)
          else
            if (zeta(n) < -zetam) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(-zetam*obu(n)/(10._r8+z0m(n))) ) )
            else if (zeta(n) < 0._r8) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n))) ) )
            else if (zeta(n) <=  1._r8) then
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n))) ) )
            else
              u10_clm(pp) = um(n) - ( ustar(n)/vkc*(log(obu(n)/(10._r8+z0m(n))) ) )
            end if
          end if
          va(pp) = um(n)
        end do
      else
        if (zldis(n)-z0m(n) .le. 10._r8) then
          u10_clm(n) = um(n)
        else
          if (zeta(n) < -zetam) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(-zetam*obu(n)/(10._r8+z0m(n))) ) )
          else if (zeta(n) < 0._r8) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n))) ) )
          else if (zeta(n) <=  1._r8) then
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(zldis(n)/(10._r8+z0m(n))) ) )
          else
            u10_clm(n) = um(n) - ( ustar(n)/vkc*(log(obu(n)/(10._r8+z0m(n))) ) )
          end if
        end if
        va(n) = um(n)
      end if

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_t_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_t_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp1(n)=vkc/log(-zetat*obu(n)/z0h(n))
      else if (zeta(n) < 0._r8) then
         temp1(n)=vkc/log(zldis(n)/z0h(n))
      else if (zeta(n) <= 1._r8) then
         temp1(n)=vkc/log(zldis(n)/z0h(n))
      else
         temp1(n)=vkc/log(obu(n)/z0h(n))
      end if

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_q_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_q_pft(n)-displa(n)
      end if
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp2(n)=vkc/log(-zetat*obu(n)/z0q(n))
      else if (zeta(n) < 0._r8) then
         temp2(n)=vkc/log(zldis(n)/z0q(n))
      else if (zeta(n) <= 1._r8) then
         temp2(n)=vkc/log(zldis(n)/z0q(n))
      else
         temp2(n)=vkc/log(obu(n)/z0q(n))
      end if

      zldis(n) = 2.0_r8 + z0h(n)
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp12m(n)=vkc/log(-zetat*obu(n)/z0h(n))
      else if (zeta(n) < 0._r8) then
         temp12m(n)=vkc/log(zldis(n)/z0h(n))
      else if (zeta(n) <= 1._r8) then
         temp12m(n)=vkc/log(zldis(n)/z0h(n))
      else
         temp12m(n)=vkc/log(obu(n)/z0h(n))
      end if

      zldis(n) = 2.0_r8 + z0q(n)
      zeta(n) = zldis(n)/obu(n)
      if (zeta(n) < -zetat) then
         temp22m(n)=vkc/log(-zetat*obu(n)/z0q(n))
      else if (zeta(n) < 0._r8) then
         temp22m(n)=vkc/log(zldis(n)/z0q(n))
      else if (zeta(n) <= 1._r8) then
         temp22m(n)=vkc/log(zldis(n)/z0q(n))
      else
         temp22m(n)=vkc/log(obu(n)/z0q(n))
      end if
      ! diagnose 10-m wind for dust model (dstmbl.F)
      ! Notes from C. Zender's dst.F:
      ! According to Bon96 p. 62, the displacement height d (here displa) is
      ! 0.0 <= d <= 0.34 m in dust source regions (i.e., regions w/o trees).
      ! Therefore d <= 0.034*z1 and may safely be neglected.
      ! Code from LSM routine SurfaceTemperature was used to obtain u10

      if (present(landunit_index)) then
        zldis(n) = forc_hgt_u_pft(pfti(n))-displa(n)
      else
        zldis(n) = forc_hgt_u_pft(n)-displa(n)
      end if 
      zeta(n) = zldis(n)/obu(n)
      if (min(zeta(n), 1._r8) < 0._r8) then
         tmp1 = (1._r8 - 16._r8*min(zeta(n),1._r8))**0.25_r8
         tmp2 = log((1._r8+tmp1*tmp1)/2._r8)
         tmp3 = log((1._r8+tmp1)/2._r8)
         fmnew = 2._r8*tmp3 + tmp2 - 2._r8*atan(tmp1) + 1.5707963_r8
      else
         fmnew = -5._r8*min(zeta(n),1._r8)
      endif
      if (iter == 1) then
         fm(n) = fmnew
      else
         fm(n) = 0.5_r8 * (fm(n)+fmnew)
      end if
      zeta10 = min(10._r8/obu(n), 1._r8)
      if (zeta(n) == 0._r8) zeta10 = 0._r8
      if (zeta10 < 0._r8) then
         tmp1 = (1.0_r8 - 16.0 * zeta10)**0.25_r8
         tmp2 = log((1.0_r8 + tmp1*tmp1)/2.0_r8)
         tmp3 = log((1.0_r8 + tmp1)/2.0_r8)
         fm10 = 2.0_r8*tmp3 + tmp2 - 2.0_r8*atan(tmp1) + 1.5707963_r8
      else                ! not stable
         fm10 = -5.0_r8 * zeta10
      end if
      if (present(landunit_index)) then
        tmp4 = log( max( 1.0_r8, forc_hgt_u_pft(pfti(n)) / 10._r8 ) )
      else
        tmp4 = log( max( 1.0_r8, forc_hgt_u_pft(n) / 10._r8 ) )
      end if
      if (present(landunit_index)) then
        do pp = pfti(n),pftf(n)
          u10(pp) = ur(n) - ustar(n)/vkc * (tmp4 - fm(n) + fm10)
          fv(pp)  = ustar(n)
        end do 
      else
        u10(n) = ur(n) - ustar(n)/vkc * (tmp4 - fm(n) + fm10)
        fv(n)  = ustar(n)
      end if
   end do

#endif

   end subroutine FrictionVelocity

!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: StabilityFunc
!
! !INTERFACE:
   real(r8) function StabilityFunc1(zeta)
!
! !DESCRIPTION:
! Stability function for rib < 0.
!
! !USES:
      use shr_const_mod, only: SHR_CONST_PI
!
! !ARGUMENTS:
      implicit none
      real(r8), intent(in) :: zeta  ! dimensionless height used in Monin-Obukhov theory
!
! !CALLED FROM:
! subroutine FrictionVelocity in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
!
!
! !LOCAL VARIABLES:
!EOP
      real(r8) :: chik, chik2
!------------------------------------------------------------------------------

      chik2 = sqrt(1._r8-16._r8*zeta)
      chik = sqrt(chik2)
      StabilityFunc1 = 2._r8*log((1._r8+chik)*0.5_r8) &
           + log((1._r8+chik2)*0.5_r8)-2._r8*atan(chik)+SHR_CONST_PI*0.5_r8

    end function StabilityFunc1

!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: StabilityFunc2
!
! !INTERFACE:
   real(r8) function StabilityFunc2(zeta)
!
! !DESCRIPTION:
! Stability function for rib < 0.
!
! !USES:
     use shr_const_mod, only: SHR_CONST_PI
!
! !ARGUMENTS:
     implicit none
     real(r8), intent(in) :: zeta  ! dimensionless height used in Monin-Obukhov theory
!
! !CALLED FROM:
! subroutine FrictionVelocity in this module
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
!
!
! !LOCAL VARIABLES:
!EOP
     real(r8) :: chik2
!------------------------------------------------------------------------------

     chik2 = sqrt(1._r8-16._r8*zeta)
     StabilityFunc2 = 2._r8*log((1._r8+chik2)*0.5_r8)

   end function StabilityFunc2

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: MoninObukIni
!
! !INTERFACE:
  subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu)
!
! !DESCRIPTION:
! Initialization of the Monin-Obukhov length.
! The scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
! Vol. 11, 2628-2644.
!
! !USES:
    use clm_varcon, only : grav
!
! !ARGUMENTS:
    implicit none
    real(r8), intent(in)  :: ur    ! wind speed at reference height [m/s]
    real(r8), intent(in)  :: thv   ! virtual potential temperature (kelvin)
    real(r8), intent(in)  :: dthv  ! diff of vir. poten. temp. between ref. height and surface
    real(r8), intent(in)  :: zldis ! reference height "minus" zero displacement heght [m]
    real(r8), intent(in)  :: z0m   ! roughness length, momentum [m]
    real(r8), intent(out) :: um    ! wind speed including the stability effect [m/s]
    real(r8), intent(out) :: obu   ! monin-obukhov length (m)
!
! !CALLED FROM:
! subroutine BareGroundFluxes in module BareGroundFluxesMod.F90
! subroutine BiogeophysicsLake in module BiogeophysicsLakeMod.F90
! subroutine CanopyFluxes in module CanopyFluxesMod.F90
!
! !REVISION HISTORY:
! 15 September 1999: Yongjiu Dai; Initial code
! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision
!
!
! !LOCAL VARIABLES:
!EOP
!
    real(r8) :: wc    ! convective velocity [m/s]
    real(r8) :: rib   ! bulk Richardson number
    real(r8) :: zeta  ! dimensionless height used in Monin-Obukhov theory
    real(r8) :: ustar ! friction velocity [m/s]
!-----------------------------------------------------------------------

    ! Initial values of u* and convective velocity

    ustar=0.06_r8
    wc=0.5_r8
    if (dthv >= 0._r8) then
       um=max(ur,0.1_r8)
    else
       um=sqrt(ur*ur+wc*wc)
    endif

    rib=grav*zldis*dthv/(thv*um*um)
#if (defined PERGRO)
    rib = 0._r8
#endif

    if (rib >= 0._r8) then      ! neutral or stable
       zeta = rib*log(zldis/z0m)/(1._r8-5._r8*min(rib,0.19_r8))
       zeta = min(2._r8,max(zeta,0.01_r8 ))
    else                     ! unstable
       zeta=rib*log(zldis/z0m)
       zeta = max(-100._r8,min(zeta,-0.01_r8 ))
    endif

    obu=zldis/zeta

  end subroutine MoninObukIni

end module FrictionVelocityMod
